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Abstract 


A multiporosity extension of classical double and triple porosity fractured rock flow 
models for slightly compressible fluids is presented. The multiporosity model is an 


adaptation of the multirate solute transport model of Haggerty and Gorelick (1995) to 


viscous flow in fractured rock reservoirs. It is a generalization of both pseudo-steady- 
state and transient interporosity flow double porosity models. The model includes a 
fracture continuum and an overlapping distribution of multiple rock matrix continua, 
whose fracture-matrix exchange coefficients are specihed through a discrete probability 
mass function. Semi-analytical cylindrically symmetric solutions to the multiporosity 
mathematical model are developed using the Laplace transform to illustrate its behav¬ 
ior. The multiporosity model presented here is conceptually simple, yet flexible enough 
to simulate common conceptualizations of double and triple porosity flow. This com¬ 
bination of generality and simplicity makes the multiporosity model a good choice for 
flow in low-permeability fractured rocks. 


1 Introduction 

The flow of slightly compressible fluids through fractured rocks is of fundamental importance 
to groundwater and hydrocarbon production and effective isolation of radioactive waste and 
CO2. In low-permeability rocks, fractures are often the primary source of bulk permeability. 
Reservoir and storage rocks often include multiple overlapping and intersecting fracture 
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sets. Discrete representations of fractures in numerical models is possible if information is 
available about fracture spacing, orientation, and aperture. Typically, spatial distributions 
of these data are not available, and fractured rocks must be approximated as porous media 

Multiple porosity types can simultaneously participate in 


Gringarten 

1982 

Sahimi 

2011) 


flow through fractured and unfractured rocks. Shales have macroscopic porosity associated 


with clastic particles and microscopic porosity associated with their organic fraction (Akkutlu 


and Fathi, 2012). The Culebra Dolomite has illustrated effects of multiple porosity types 


(fractures, vugs, and intergranular porosity) during flow and solute transport at both the 
held and laboratory scale (McKenna et ah, 2001 Malama et ah, 2013). The effect of multiple 
types and scales of porosities can be masked by or mistaken for the effects of multiple types 


and scales of heterogeneity (Altman et ah, 2002). 


Dual or multiple porosity models are common simplihed conceptualiztaions of this com¬ 
plex reality of interacting porosity types and the spatially heterogeneous nature of rocks. The 
treatment of fractured rock as a system of interacting and overlapping continua has been 


1988 

Chen| 

1989; 

Da Prat 

1990 

; or Bourdet 

2002) 


2002) since its hrst introduction by Baren- 


blatt and Zheltov (1960). The conceptualization is based on a low-permeability but high 


storage-capacity rock matrix drained via natural and man-made high-permeability but low 
storage-capacity fractures. 

Slightly compressible huids include water, brine, oil, and other liquids for which density 
can be assumed constant, allowing a formulation in terms of head (rather than pressure). 
Gases can be treated as a slightly compressible huids after a transformation to create a 
pseudo-pressure and pseudo-time, which accounts for pressure-dependence of huid proper¬ 
ties (Friedmann, 1958 A1 Hussainy et al., 1966). Slightly compressible how models can also 
approximate more complex non-linear how (i.e., multiphase how, time-variable permeabil¬ 
ity, or gas desorption) by transforming observations into one of several possible integrated 
pseudo variables (see review by Clarkson, 2013). Flow prediction for slightly compressible 
huids in fractured low-permeability rocks at a macroscopic (wellbore or reservoir) scale is 
of great importance to applications in groundwater supply, hydrocarbon production, and 
underground sequestration of nuclear waste or CO 2 . 

Fractured reservoirs exhibit pressure drawdown due to production at a specihed pressure 
or production at a specihed howrate (e.g., during recovery or shut-in), in a manner char¬ 
acteristic of multiple interacting porosities (e.g., Crawford et ah, 1976; Gringarten, 1982 


Moench, 1984). Fractures provide high-permeability pathways, often orders of magnitude 


more permeable than the unfractured rock itself, but comprise only a small portion of the 
rock volume. Natural fracture porosity is often less than 0.1%, but the fracture network is a 
much more efficient huid conductor per unit porosity than the matrix (Streltsova, 1988). Sev¬ 
eral well-known double porosity conceptualizations for “uniformly” or “naturally” fractured 


Barenblatt et al.| 


Warren and Root 

1963 


of fractures). These solutions represent the how domain with overlapping fracture and ma¬ 
trix continua. These models assume the domain of interest is large enough and the fracture 
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density is uniform enough to treat both fractures and matrix as continua. Fractures are the 
main source of permeability and spatial connectivity, while the matrix is the main source of 
storage. Triple porosity solutions with a single fracture continuum and two matrix continua 
(only connected through the fracture porosity) are a logical extension of dual porosity (e.g., 


Clossman 1975 Liu 1981 |A1 Ahmadi and Wattenbarger, 2011 and Tivayanonda et ah 


2012 ). 


Beginning from the system proposed by Barenblatt and Zheltov (1960), but equally al¬ 


lowing advection through fractures and matrix results in the dual permeability model, which 
requires a coupled set of governing equations and boundary conditions. These were hrst 


solved analytically by Chen and Jiang (1980) in terms of quasi-Bessel functions, but this so¬ 


lution has not seen wide use. Although the dual-permability approach is more physically re¬ 
alistic than the simplihed dual-porosity approach, when matrix permeability is much smaller 


than the fractures, dual porosity is an adequate approximation (Chen, 1989). The multiple 


interacting continua (MINC) approach taken by TOUGH2 is a numerical approach capable 


of simulating both dual-porosity and dual-permeability systems (jPruess and Narasimhan 
1985 Pruess et al.| 1999) 


Many analytical solutions for variations on double and triple porosity conceptual models 


have been developed (Moench, 1984 Chen, 1989 Da Prat, 1990), with different conhgu- 


rations and relationships between fracture and matrix continua, but a generalized solution 
is needed. We present multiporosity (the name coming from multirate and dual-porosity) 
as a generalization and extension of the Warren and Root ( |1963 ) solution to any number 
of matrix continua interacting with a fracture continuum. Multiporosity is an adaptation 
of the multirate solute transport theory to viscous flow, for computing pressure-driven flow 
through low-permeability rocks with heterogeneity and multiple porosities. The multirate 


(i.e., multiple reaction rates) conceptual model for solute transport (Haggerty and Gorelick 


1995, 1998) has been successfully used to simulate diffusion of solutes from fast-flowing frac¬ 


tures into a distribution of diffusion-dominated matrix block sizes (Haggerty et ah, 2000 


2001 

McKenna et ah 

2001 

Malama et al. 

2013 


2 Multiporosity Model 

It is useful to conceptualize slightly compressible flow in uniformly fractured domains as dif¬ 
fusion in a porous medium. Some deviations from this ideal behavior can be accommodated 


through use of pseudo-time and -pressure (Clarkson, 2013) 


2.1 Conceptualization 

Two primary flow conceptualizations used in slightly compressible dual-porosity systems are: 


1. The pseudo-steady-state Warren and Root (1963) (WR) interporosity flow concep¬ 


tualization (similar to Barenblatt et al(| 1960) assumes flow from matrix blocks is 


proportional to the pressure difference between the fracture and the average matrix 
pressures. 
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2. The Kazemi (1969) (KZ) transient interporosity flow conceptualization allows transient 
diffusion from the fracture to the matrix, and couples flow in the fracture and matrix 
domains through a source term proportional to flux in the fracture flow governing 
equation. 

The WR matrix flow conceptualization is much simpler than the KZ conceptualization, 
but leads to analytical flow solutions more readily than the KZ approach (e.g., see review by 


Chen 


The KZ approach has been solved using finite differences (Kazemi, 1969), finite 


1985, or Ozkan et ah, 1987) and approximate (Walton and McLennan 


volumes ([Pruess et ah 1999), and various analytical (e.g., Serra et ah 


1983 


Chen et al. 


2013) approaches. 


The WR solution produces physically unrealistic pressure transient solutions during the 
transition from early fracture-dominated flow to later matrix-dominated flow (Gringarten 


1982 Moench, 1984). It produces a nearly flat transition between early fracture flow and 


late-time matrix flow on a semi-log plot, equivalent to a vanishing log time derivative - 
d/d{\at). Moench (1984) showed the WR model could be improved by adding a fracture skin. 


which delayed communication between the fracture and matrix. The WR approach assumes 
pseudo-steady-state flow between fracture and matrix, which may not occur physically until 
late time in low-permeability rocks. The KZ model produces a half-slope transition between 
the early and late-time flow, which is more in agreement with typical held observations 


(Moench, 1984). 


The multiporosity approach is a generalization of the WR double-porosity model and the 


pseudo-steady-state triple-porosity model of Clossman (1975). Triple-porosity models can 


represent two systems of natural preexisting micro- and macro-fractures (Al Ahmadi and 


Wattenbarger, 2011). The multiporosity model is presented here as a spatial distribution of 


natural fractures and matrix materials (i.e., matrix heterogeneity). 

When moving between models with different numbers of porosities (e.g., double- or triple¬ 
porosity compared to single-porosity models), a model with more parameters is typically 
more flexible and able to fit a wider range of observed behaviors, but more free parameters 
must be estimated. Data collected from typical low-permeability wells are inadequate to 
uniquely constrain models with many estimable parameters. To constrain the number of 
free parameters in the multiporosity model, the interporosity flow parameters can be speci¬ 
fied with a distribution function (e.g., Ranjbar et al. (2012[)). Haggerty and Gorelick (1995) 


derived a distribution comprised of an infinite series of porosities. Their distribution is equiv 
alent to a distribution of pseudo-steady-state WR matrix porosities that behave in total as 
a transient KZ-type system. The multiporosity model can be made equivalent to diffusion 
into a slab (i.e., the KZ conceptualization), cylinder, or sphere. During parameter estima¬ 
tion, the multiporosity model can be matched to data either with flexible but parsimonious 
property distributions (e.g., lognormal or beta), or with an arbitrary number of individual 
porosities and their associated parameters. The multiporosity distribution can be flexible, 
but with certain distributions it can also be shown to be a generalization of two well-known 
physically based end members. 

Porosity distributions have been utilized in some dual-porosity solutions based upon dif¬ 
ferent block-size distribution function 


McGuinness 

l').S(i 

Ghen 


and 

Ranjbar 
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et al., 2012), which is similar to the approach taken here. Through parameter estimation, 


matrix block distributions can provide a fracture distribution, which is much more flexible 
than typically double or triple porosity models that assume more uniformly spaced fracture 
distributions. Existing block-size distribution solutions have focused primarily on the deriva¬ 
tion of shape factors. We show how certain distributions of porosities lead to the well-known 
special cases of WR and KZ double-porosity flow. 

A multiporosity conceptual model for flow has two possible meanings. Multiporosity can 
represent a single matrix continuum with the inter-porosity exchange coefficient between 
the fracture and matrix continua treated as heterogeneous (i.e., a random variable). This 
is a simplification of small-scale spatial heterogeneity (e.g., Haggerty and Gorelick, 1995 


Haggerty and Gorelick, 1998 or Akkutlu and Fathi, 2012). Alternatively, it can represent 


a set of multiple physically-distinct matrix continua with the porosity and permeability of 
each as random variables. Our discussion takes the former approach, we believe it to be the 
most physically realistic. 


2.2 Mathematical Model 


We present the multiporosity flow model, which is a logical extension of the WR conceptu¬ 


alization (Warren and Root, 1963), to an arbitrary number of pseudo-steady-state matrix 


domains. The model is an adaptation of the multirate advection-dispersion solute transport 
solution (Haggerty and Gorelick, 1995, 1998) to viscous flow. It conceptualizes random poros¬ 
ity spatial variability in rock matrix as a random distribution of uniform matrix porosities 
communicating with the primary fracture porosity. 

The governing equation for pressure head drawdown Ap(x, f) = po(x) — p(x, f) (po is 
initial pressure head) in the fracture continuum is 




d (Ap/) 

dt 


N 


j=0 


d{Apj) 

dt 




(Apf ), 


( 1 ) 


where t is time, index j and subscript / denote quantities related to the jth matrix and frac¬ 
ture continua, the sum is across N matrix porosities {N may be infinite), 0 is dimensionless 
porosity, c is a compressibility or storage coefficient, k is permeability, p is fluid viscosity, 
and X is a dimensionless probability mass function (PMF - the discrete form of a probabil¬ 
ity density function) of interporosity exchange coefficients. See Table for a summary of 
physical quantities and their units. 

In this multiporosity conceptualization, matrix properties and dependent variables are 
implicitly discrete functions of the distribution of matrix continua index. The properties 
controlling matrix-fracture fluid exchange across a potentially infinite number of matrix- 
fracture interfaces can be considered a random variable, but the resulting governing equations 
are deterministic because only their sum or bulk behavior appears in ([^. 

Flow in each of the matrix domains is generally governed by the diffusion equation, viz. 

= WyAp,) j = (2) 
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the governing equations will be non-dimensionalized and solved using the Laplace transform. 

The governing multiporosity fracture flow equation Q can be non-dimensionalized by 
dividing through by the total formation compressible storage (j)c = (jyfCf + 4’jCj-, and a 

characteristic pressure Pc, to produce the dimensionless fracture flow equation 


Uf 


dtn 


j=i 




(3) 


where 'ijji = Ap^/Pc {i G {/, j}) is dimensionless pressure change, Ui = (i)e,ciil {cj)c) is the 
fractional storage of an individual continuum (0 < < 1 and ~ 1)) VI/Ll = 

is the dimensionless Laplacian, Lc is a characteristic length (often dehned as the pumping 
well radius, r^, to = t/Tc, and Tc = L‘/.p(f)c/kf is a characteristic time. Alternatively, Ui 
can be related to the volume-weighted storage coefficient commonly used in h ydrogeology, 
= SsjVi/ (^SgfVf + Yj ^sjVj^, where Id is a fraction of the total volume (Gringarten 


1982 Moench, 1984|). Table defines dimensionless quantities for the current problem 


The matrix flow equation (|^ can analogously be non-dimensionalized into 


= K,jV 


dt 


D 


ni’j 




(4) 


where Kj = kj/kf is the jth matrix-to-fracture permeability ratio. 

To simplify from the transient to the pseudo-steady-state interporosity flow conceptu¬ 
alization we use average matrix pressure by integrating the matrix governing equation (|^ 
across the matrix domain. We assume matrix blocks are comprised of one-dimensional 
slabs (i.e., rectangular radially symmetric blocks with flow perpendicular to fractures, see 
Figure [^. The layered idealization in Figure results in an equivalent solution to the 


dual porosity conceptualization when continua exist at each physical location (Streltsova 


1988). We choose slabs for simplicity; other possible matrix block geometries (e.g., spheres 


or cylinders) do not result in markedly different predicted results for the double porosity 
conceptualization (e.g., Gringarten, 1982; Moghadam et ah 1 120101). Integrating (|4l) across a 


one-dimensional slab results in 

d{^j) 

Kj 



d^j 

1 

dto 

Lnooj 

dyo 

yD=Ln 

dyo 

yD=o_ 


(5) 


where Ljj = L/Lc is the dimensionless half-distance between evenly spaced fractures, pd = 
y/Lc is the dimensionless matrix space coordinate perpendicular to the fracture (0 < < 

Lb), and 

1 

Jo 

is the spatially averaged dimensionless change in matrix pressure. The first term on the 
right-hand-side of ([^ is flux at the matrix-fracture interface, while the second term is flux 
at the center of the block. The latter vanishes identically by symmetry. 
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The two no-flow boundaries perpendicular to the borehole in Figure are symmetry 
boundary conditions, which deflne the unit cell represented by the model given here. In the 
case of a systematically fractured horizontal well one plane is in the midplane of the fracture 
and the other halfway between equally spaced fractures. The cylindrical boundary parallel 
to the borehole represents the edge of the reservoir - potentially at infinite radial distance. 
At this far edge of the domain we implement a general linear Type-III boundary condition. 

Continuing with the WR pseudo-steady-state interporosity flow assumption, flux from 
each matrix continuum to the fracture continuum is proportional to the difference between 
the dimensionless fracture and average matrix pressure changes. 


dvD 


yD=LD 



[V'z - ii’j)] 


3 = 1, 


,iV, 


( 6 ) 


where Cj is a dimensionless constant of proportionality. Substituting this expression for flux 
0 into the integrated matrix flow equation ([^ leads to 


dto 


LjyUj 


[^/ - (^j)] 




(7) 


Comparing the j = 1 matrix porosity from ([t]) t o the standard WR solution for matrix flow 
in a double-porosity system (Warren and Root, 1963, Eqn. 11), results in the equivalence 


arlkm 


eikiL 


{1 — u})kf L'^Uikf' 


( 8 ) 


which suggests ei = where a is WR’s shape parameter [L“^], k^ = ki is the WR 

matrix permeability, oj = ojf, oji = 1 — oj, and WR used Lc = as a characteristic length. 
We choose the group 


to characterize the flow problem across the distribution of matrix continua (0 < Uj < oo). 
Taking the Laplace transform of ([^ results in 

{^j) s = Uj [ijf - {^j)] j = l,...,N (9) 


where s is the dimensionless Laplace transform parameter and an overbar indicates a trans¬ 
formed dependent variable, i.e., / = flin) dtf). The averaged change in dimension¬ 

less matrix pressure due to changes in the fracture pressure is 


(^i> 


S + Uj 


This can be substituted into the Laplace-transformed form of (|^ after similarly integrating 
(|^ across the matrix blocks (equivalent to replacing tjjj with {i/jj)), resulting in 


N 

SCJfijf + SCJf 

i=i 


/^jXj 


Ui 


S - 1 - Uj 


V’/ = V^V’/, 


( 10 ) 
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where I3j = ujj/ojf a. dimensionless matrix/fracture storage capacity ratio. 
Equation (10) can be further simplihed into 


{l + g) = 0 , 


where 


9 = 


N 

E 


XjUj 

S + M, ’ 


( 11 ) 


( 12 ) 


is the matrix memory kernel (|Haggerty and Gorelick, 1995) and Xj = Xjf^j is a scaled PMF. 


2.3 Distributions 

To compute a solution to ( 0 . the matrix memory kernel must be specihed. Any valid 
PMF can be specified for Xj, but we present three special cases. The multiporsity solution 
simplihes to the dual-porosity WR solution through 


5 'wr — 


- 1) A 

\UJ / 1—LJ 


X 

1 — LU 


(13) 


which is Ml = A/(l — ca) = X/ui = xi = 1, and N = 1. Here A = is WR’s 

dimensionless interporosity flow parameter. 

Similarly, the multiporosity solution is equivalent to the pseudo-steady-state triple-porosity 
solution of Clossman (1975) or Liu (1981| using (12), 


7i(l -7) 


Xi = 


7 


and N = 2, where 7 is the bulk volume fraction of hssures, an d 7 ( 1 , 2 } are the volume fraction 
of good and poor rock in the matrix (Odeh, 1965 Clossman, 1975). 


Haggerty and Gorelick (1995, Table 1) presented inhnite discrete distributions which 


make the overall multiporosity system (comprised of a sum of WR pseudo-steady-state con- 
tinua) behave like transient interporosity flow (i.e., the KZ model). They presented series of 
coefficients for a matrix memory kernel, which are mathematically equivalent to diffusion into 
matrix blocks of different geometry through analogies between the series solution and ana¬ 


lytical solutions (e.g., similar to those for dual porosity by de Swaan O. (1976) and Najurieta 


(1980)). Diffusion into a slab is equivalent to an inhnite distribution of pseudo-steady-state 
matrix domains given by 


{2i - lyTT^eiKi . 

Ui = -^- ? = 0 , 1 , 


and 


Xi = 


Auji 


8ui 


{2i — lyufTr'^ 




(14) 


(15) 









































Haggerty and Gorelick (1995) presented similar series of coefficients representing diffusion 
into cylinders or spheres (as did de Swaan O., 1976). For their problem, they found truncating 


the series at iV = 100 resulted in pseudo-steady-state solutions within 1% of matrix diffusion, 
with errors conhned to early time. The series are quick to compute, and using 1000 or more 
terms is not computationally expensive. 

Alternative distributions of matrix continua and properties can be specified by other 
widely used probability distributions, such as exponential, normal, lognormal, or linear (e.g.. 


Ranjbar et ah 

2012 

Malama et ah, 

2013 


3 Solutions 


The solution to (11) in cylindrical coordinates can be found directly in Laplace space, analo¬ 


gous to solutions in the literature for traditional dual-porosity problems (Warren and Root 


1963 Kazemi, 1969| Mavor and Cinco-Ley, 1979) for either Type I (specihed down-hole 


pressure) or Type II (specified flowrate) wellbore boundary conditions. 


3.1 Flow problem of interest 

The following assumptions and boundary conditions are used to develop Laplace-space an¬ 
alytical solutions to the governing equations derived in the previous section: 

1 . single-phase slightly compressible flow (i.e., water, brine, oil, or gas treated using 
appropriate pseudo-variable methods), 

2. a specified flowrate Q(t) (Type II) or pressure change Ap(t) (Type I) at the well 
completion, 

3. the completion only intersects or interacts with fractures, with the matrix connected 
to the completion through fractures, 

4. a symmetry no-flow boundary conditions parallel to the fracture midway between two 
fractures, and 

5. a Type-Ill boundary condition at far edge of the domain, r = R (Figure [^. The Type- 
III boundary condition can either represent no-flow (bounded reservoirs), specified head 
(a circular island or laboratory sand tank with a specified head condition), or a linear 
combination of the two. 


The pressure solution is developed for drawdown, ip = {p — po) /Pc (change from an initial 
state), therefore only the homogeneous initial condition is considered. 

The Laplace-transformed governing flow equation (0 in cylindrical coordinates with 
radial symmetry is the modihed Helmholtz equation for radial coordinates. 


d'^ipf I dipf 


0 , 


(16) 
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where r^) = r/L^ and rf = sujf (1 + g) is a purely imaginary wave number, also called the 
fracture function (A1 Ahmadi and Wattenbarger, 2011). We pick Lc = as a. characteristic 


length for radial flow to a well. 


3.2 Approach 


Analogous to Mavor and Cinco-Ley (1979), we compute time-domain numerical values from 
solutions to (16) using a numerical inverse Laplace transform algorithm (de Hoog et al. 


1982). Several feasible alternative numerical inversion approaches exist (Kuhlman, 2013), 


some which only require real values of s; the de Hoog et al.| (1982) approach requires complex 


s. 


We present solutions for both hnite and inhnite domains as well as for specihed bottom- 
hole pressure or flowrate using modihed Bessel functions - the radial eigenfunctions of the 
modihed Helmholtz equation in cylindrical coordinates. 

The multiporosity conceptualization may be appropriate for flow in fractured rock in 
which hydraulic fracturing is being performed, but the simple homogeneous cylindrical geom¬ 
etry of the radial solution presented here cannot adequately represent pumping discrete frac¬ 
tures under typical held conditions (GringartenH 1982 ,[crarkson 2013). The governing equa¬ 


tions presented here could be solved for discrete linear fractures using analytical solutions for 


single fractures (e.g.. 

Gringarten et al. 

1974 

) or using analytic element superposition-based 

combinations of line element solutions 

'Bakker and Kuhlman 

2011 ; 

Biryukov and Kuchuk 


2012). These approaches can better represent the geometry encountered during hydraulic 


fracturing. Line sinks and sources, or narrow ellipses and polygons of high permeability can 
be used to represent discrete fractures added to a uniformly fractured multiporosity rock. 


3.3 Infinite Domain Solutions 


For an inhnite radial domain we assume a homogeneous far-held Type-I boundary condition, 
lim^^-^oo 'ipfi'f'o) = 0, while the wellbore boundary condition is specihed in two diherent ways 
as follows. 

For specihed down-hole howrate, we choose Pc = fiBq/{27rhfkf), where Q(t) = qft is 
a volume howrate at the surface, g is a constant characteristic howrate, hj is the portion 
of the well completion open to fractures, and i? is a dimensionless formation volume factor 
correcting for the diherence between down-hole and surface howrates due to compressibility 
(in low-pressure or incompressible how systems B = 1). This choice of Pc results in the 
wellbore boundary condition 



dro 


ro=l 



where ft is the Laplace-transformed temporal behavior of the boundary condition (i.e., tem¬ 
poral huctuations from q). 

Using ft, we consid er general time behavior using the Laplace-space form of Duhamel’s 
theorem (Ozi§ik, 1993, Chap. 5). For example, a step function on at to = too is ft = 
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(with the typical case of ft = l/s for t^o = 0), and a pulse function during too < is 

ff = / s. The temporal behavior of an arbitrary flowrate, Q{t), or down-hole 

pressure, Ap(t), is readily approximated using piecewise constant or linear behavior (e.g.. 


Streltsova, 1988 or Mishra et ah, 2013). For nearly constant or monotonically declining 


specihed pressures and flowrates, the numerical inverse Laplace transform algorithm con¬ 
verges quickly (e.g., 21 Laplace-space function evaluations for each time is typical). Many 
more Laplace-space function evaluations are needed to ensure accuracy when using periodic 


or rapidly fluctuating rates (Kuhlman, 2013). For this reason, we recommend htting ft to 


follow general observed trends, rather than every observed fluctuation. 

The Laplace-space solution to (16) for in the completion {r^ = 1), given this boundary 
condition at the wellbore, is 




(17) 


where Kn{x) is the second-kind modihed Bessel function of order n (DLMF, §10). 

We do not consider wellbore storage explicitly. Mavor and Cinco-Ley ( 1979| ) presented 
similar solutions including both wellbore storage and skin effects; these solutions could be 
used with the multiporosity conceptualization to approximately include wellbore storage 
and skin effects. To rigorously include wellbore storage effects would additionally require 
including terms related to the time derivative of ft- 

For specihed downhole pressure, we chose Pdfo) = Poft- The Laplace-space solution to 


(16) for howrate in the completion, given the Type-I boundary condition at the wellbore, is 




(18) 


3.4 Finite Domain Solutions 


For a hnite cylindrical domain centered on the well completion, with a homogeneous Type- 


Ill boundary condition at r = R {td = Rd)-, the Laplace-space solution to (16) is given 


generally by Carslaw and Jaeger ( 1959| §13.4). The specihc Laplace-space solution to (16) 
for m the completion, given the specihed howrate boundary condition at the wellbore 
and the Type-III boundary condition + HApf = 0 at r = R (nondimensionalized in 

Laplace space to -|- = 0, where Hd = HpLc/kf), is 


Tiq,R) ^ ft foiv)^ + f^oiv)C 

f h - Ki(77)C 


(19) 


where 


^ = riKiij^RD) - HoRoiTlRD), 
C = rj\i{'qRD) + HDfoi'qRD), 


H is the Type-III boundary surface conductivity, and h 
function of order n (DLMF, §10). 


(x) is the hrst-kind modihed Bessel 
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The Laplace-space solution to (16) for flowrate, given the specihed pressure boundary- 


condition at the wellbore and a homogeneous Type-III boundary condition at the outer 
boundary, is 


^ -irh/Drift 


h(v)( - KiMC 


( 20 ) 


In the limiting case where Hd = 0, the boundary condition sX r d = Ro becomes a no-flow 
condition. When Hd ^ oo, the outer boundary condition becomes a specihed-pressure 
condition. In the more general case, Hd represents “leakiness” of the outer boundary; a 
second reservoir beyond td = Rd is providing a non-zero flux into the domain, proportional 
to the change in pressure at the boundary. 


4 Model Behavior 


We present examples of some typical behavior of the multiporosity model. First we illustrate 
how the series of matrix porosities (Equations [l4| and [l5| can approximate both WR and KZ 
flow as end members for the inhnite domain solution. Second, we illustrate model predictions 
for a range of different matrix and fracture properties. Finally, we illustrate the nature of the 
different boundary conditions {H G {0, if, cxd}) at the radial extent of the domain td = Rd- 

For constant matrix properties, the multiporosity model can reproduce both pseudo- 
steady-state (WR) and transient dual-porosity (KZ) flow, as well as a range of intermediate 
behaviors. Figure shows the behavior of the inhnite domain solution 0 for Uj = u, 
/3j = 13, and different values of N. Using constant formation properties for an inhnite number 
of pseudo-steady-state matrix porosities results in an equivalent transient KZ solution with 
similar properties. The hgure shows the increase in dimensionless pressure drawdown due 
to production at a constant rate. The line types indicate the model, while the line colors 
indicate the matrix/fracture permeability ratio k. The right plot in Figure shows the 
slope with respect to ln(t), illustrating how the KZ model (solid line) has an intermediate 
slope which is 1/2 of its slope at early or late time. The WR model has essentially zero 
slope at intermediate time. As more terms are added to the multiporosity model, the slope 
increases from zero to 1/2 the late time slope, beginning from the later-time portion of the 
intermediate time portion and moving back towards earlier time. 

Figure shows the inhnite domain solution for howrate, given a specihed bottomhole 
pressure (18). This hgure shows the decay in howrate, due to production at a specihed 


bottomhole pressure, illustrating the variation between the WR and KZ endmembers. 

Figure shows multiporosity solutions with diherent matrix properties for each matrix 
continuum, rather than uniform properties, as shown in Figures and These solutions are 
analogous to classical double- and triple-porosity models, where the two matrix permeabili¬ 
ties represent “good” or “bad” quality rock (Clossman, 1975), or possibly heterogeneity in the 
matrix system (A1 Ahmadi and Wattenbarger, 2011). Table shows the values of Uj and Xj 
used in (12). These values are derived from ojf = 0.001 and Kj = {0.5, 0.05,5x10“^, 5x 10“®}. 


Similar to Figure for constant parameters across the diherent porosities, the hat transition 
regions beyond the hrst transition disappear with increasing number of matrix continua for 
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this selection of matrix properties. The double-porosity model has one transition zone, the 
triple-porosity model has two, and the quad-porosity has three transition periods (a single 
porosity model has no transitions). 

Figure shows the effects of boundary conditions at the radial extent of the domain. 
The solid curves represent the solution for if = 0 (a no-flow boundary condition). For 
the smaller domain (red lines), the effects of the boundary are seen before the effects of 
the matrix porosity. In larger domains (green and blue), the effects of the boundary aren’t 
observed until later time. The black curve shows the solution for the inhnite domain. The 
dotted lines represent the solution for if —)• cx) (implemented as 10®), which is effectively a 
no-drawdown condition at the edge of the domain. 

Type I and III boundary conditions are less applicable to bounded geologic reservoirs, 
but these solutions may be useful in other applications, where a high-permeability reservoir 
surrounds a double- or multiple-porosity medium. The linear relationship between flux and 
change in pressure between the domain and boundary condition at Rd is analogous to the 
pseudo-steady-state assumption for flux between the fracture and matrix. 


5 Summary 

The multiporosity model is shown to be a generalization of pseudo-steady-state double and 


triple porosity solutions of Warren and Root (1963), Clossman (1975), and Liu (1983), as 


well as the transient double porosity solution of Kazemi (1969), given different distributions 
of matrix properties. Previous solutions for interporosity flow with variable block sizes 


'Ranjbar et al., 

2012 

) have not included the connection to the solution of 

Kazemi 

(1969), 


and have focused on the block-size distribution properties. 

The cylindrically symmetric solutions presented here for continuously fractured and dis¬ 
cretely fractured systems show how the multiporosity model can be presented as a unifying 
approach to solve multiple flow regimes using a single diffusion-based model, when the same 
properties are used across the different matrix continua. When different properties are 
used in each matrix porosity, the solution can produce results based upon the interaction 
of multiple physical reservoirs (i.e., macro- and micro-fractures). The radially symmetry 
Laplace-domain solutions presented are given for specihed downhole pressure and flowrate, 
with hnite or inhnite domains, and a general Type-III boundary condition at the far extent 
of the domain. 

Although there have already been many approximate and analytical solutions proposed 


for how in fractured rock (Gringarten, 1982; Chen, 1989), we propose the multiporosity 


solution as a simple yet unifying approach to include the ehects of multiple interacting 


reservoirs. The conceptual solute transport approach of Haggerty and Gorelick (1994, 1995 


1998) has been adapted to the dihusion of pressure between fracture and matrix reservoirs. 
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Figure 1: Equivalent layered representation of multiporosity problem assuming one¬ 
dimensional slab-type matrix geometry. Matrix is white cylindrical volume on either side of 
fractures. Additional fracture-matrix pairs illustrated in gray on either side of primary pair. 



Figure 2: Multiporosity flow solution for constant downhole flowrate {ft = 1/s) inhnite 
domain, showing approximation to WR {N = 1), and KZ {N —)■ cx)) solutions. 
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Figure 3: Multiporosity flow solution for constant downhole pressure {ft = 1/s) infinite 
domain, showing approximation to WR {N = 1), and KZ (A^ —)• oo) solutions. 



Figure 4: Multiporosity flow solution for constant flowrate {ft = 1/s) hnite domain (Type 
II, Hd = 0), multiple porosity solutions with different matrix properties {kj and Uj) in each 
matrix continuum. 
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Figure 5: Multiporosity flow solution (KZ, N —>■ oo) for specified flowrate (left) and downhole 
pressure (right) for a range of domain sizes and Hjy (Type I: Ho = 10®; Type II: Ho = 0). 


Table 1: Parameters used in multiporosity model shown in Figure [4| 


Model 

N 

Uj 

Uj 

Xj 

Double 

1 

{0.999} 

{0.5} 

{999} 

Triple 

2 

{0.005, 0.994} 

{100, 0.005} 

{5, 994} 

Quad 

3 

{0.005, 0.05, 0.944} 

{100, 0.1, 5x10-5} 

{5, 50, 944} 

Penta 

4 

{0.005, 0.05, 0.5, 0.444} 

{100, 0.1, 10-^ 10-Q 

{5, 50, 500, 444} 
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Table 2: Fundamental quantities 


B 

formation volume factor (formation volume/surface volume) 

- 

c 

total system compressibility 

L-i 


fracture and jth matrix domain compressibilities 

L-i 

hf 

length of well completion open to fractures 

L 

H 

Type-III boundary condition surface conductivity at R 

T-i 

Pirn 

fracture and jth matrix domain pressure heads 

L 

Po 

initial pressure head 

L 


fracture and jth matrix domain permeabilities 

L2 

L 

half width of matrix domain blocks 

L 

Lc 

characteristic length 

L 

Pc 

characteristic pressure head 

L 

Q 

time-variable volumetric flowrate at surface 

L3. t-i 

q 

constant volumetric flowrate at surface 

L3. t-i 


wellbore radius 

L 

R 

domain radius (when hnite) 

L 

s 

Laplace transform parameter 

- 

t 

time since beginning of testing or production 

T 

Tc 

characteristic time 

T 

0 

total system porosity 

- 


fracture and jth matrix domain porosities 

- 

P 

fluid viscosity 

L-T 

X 

interporosity flow parameter probability mass function 

- 


Table 3: Derived dimensionless quantities 


hfD 

/Pc 

dimensionless fracture completion length 

Hd 

= HfiLjkf 

dimensionless boundary conductivity 

Ld 

= l/l. 

dimensionless matrix domain block length 

td 

= r/Lc 

dimensionless fracture coordinate (perpendicular to wellbore) 

to 

= t/T, 

dimensionless time 

Rd 

= R/Lc 

dimensionless domain radius 

Vd 

= y/Lc 

dimensionless coordinate parallel to wellbore 

0 . 

= 

storage capacity ratio for matrix domain j 

V 

= (1 + y) 

Helmholtz wave number for multiporosity flow 

Kj 

= kj/kf 

matrix fracture permeability ratio for matrix domain j 


~ 0{/J}'3{/,i}/ (0c) 

fracture and jth matrix domain storage ratio 
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